Analisa Gempa Indonesia 2008

Latar Belakang

Sudah bukan lagi menjadi rahasia umum bahwa Indonesia berada pada rangkaian ring of fire. Ring of fire atau yang biasa disebut, Circum-Pacific Belt, adalah rangkaian gunung berapi aktif sepanjang 40.000 km yang membentang di samudera pasifik. Berdasarkan fakta tersebut, bukan hanya Indonesia memiliki jumlah gunung api aktif, tetapi juga menunjukkan bahwa Indonesia menjadi salah satu negara yang rawan terjadi gempa bumi.

Melalui analisa ini saya ingin melihat pola dan kecenderungan peristiwa gempa bumi yang terjadi di Indonesia sepanjang tahun, dari 2008 hingga Oktober 2022. Data saya dapatkan dari situs Kaggle, yang di scraped dan di atur oleh BMKG (Badan Meteorologi, Klimatologi dan Geofisika) link.

Saya menggunakan R dan RStudio dalam proses analisa ini karena melihat data yang cukup besar, dan memikirkan efisiensi untuk mendapatkan hasil yang cepat.

Memuat Libraries

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0      ✔ purrr   0.3.5 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(dplyr)
library(lubridate)

## 
## Attaching package: 'lubridate'
## 
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(ggpubr)
library(leaflet)

Mengimpor dataset

df <- read_csv("katalog_gempa.csv")

## Rows: 89218 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr   (1): remark
## dbl  (10): lat, lon, depth, mag, strike1, dip1, rake1, strike2, dip2, rake2
## date  (1): tgl
## time  (1): ot
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(df)

## # A tibble: 6 × 13
##   tgl        ot         lat   lon depth   mag remark strike1  dip1 rake1 strike2
##   <date>     <time>   <dbl> <dbl> <dbl> <dbl> <chr>    <dbl> <dbl> <dbl>   <dbl>
## 1 2008-11-01 21:02:43 -9.18  119.    10   4.9 Sumba…      NA    NA    NA      NA
## 2 2008-11-01 20:58:50 -6.55  130.    10   4.6 Banda…      NA    NA    NA      NA
## 3 2008-11-01 17:43:12 -7.01  107.   121   3.7 Java …      NA    NA    NA      NA
## 4 2008-11-01 16:24:14 -3.3   128.    10   3.2 Seram…      NA    NA    NA      NA
## 5 2008-11-01 16:20:37 -6.41  130.    70   4.3 Banda…      NA    NA    NA      NA
## 6 2008-11-01 14:47:00 -7.37  105.    18   3.3 Java …      NA    NA    NA      NA
## # … with 2 more variables: dip2 <dbl>, rake2 <dbl>

str(df)

## spc_tbl_ [89,218 × 13] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ tgl    : Date[1:89218], format: "2008-11-01" "2008-11-01" ...
##  $ ot     : 'hms' num [1:89218] 21:02:43 20:58:50 17:43:12 16:24:14 ...
##   ..- attr(*, "units")= chr "secs"
##  $ lat    : num [1:89218] -9.18 -6.55 -7.01 -3.3 -6.41 -7.37 0.1 -7.07 -3.32 -4.43 ...
##  $ lon    : num [1:89218] 119 130 107 128 130 ...
##  $ depth  : num [1:89218] 10 10 121 10 70 18 12 135 10 10 ...
##  $ mag    : num [1:89218] 4.9 4.6 3.7 3.2 4.3 3.3 4.7 4.8 2.3 3.2 ...
##  $ remark : chr [1:89218] "Sumba Region - Indonesia" "Banda Sea" "Java - Indonesia" "Seram - Indonesia" ...
##  $ strike1: num [1:89218] NA NA NA NA NA NA NA NA NA NA ...
##  $ dip1   : num [1:89218] NA NA NA NA NA NA NA NA NA NA ...
##  $ rake1  : num [1:89218] NA NA NA NA NA NA NA NA NA NA ...
##  $ strike2: num [1:89218] NA NA NA NA NA NA NA NA NA NA ...
##  $ dip2   : num [1:89218] NA NA NA NA NA NA NA NA NA NA ...
##  $ rake2  : num [1:89218] NA NA NA NA NA NA NA NA NA NA ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   tgl = col_date(format = ""),
##   ..   ot = col_time(format = ""),
##   ..   lat = col_double(),
##   ..   lon = col_double(),
##   ..   depth = col_double(),
##   ..   mag = col_double(),
##   ..   remark = col_character(),
##   ..   strike1 = col_double(),
##   ..   dip1 = col_double(),
##   ..   rake1 = col_double(),
##   ..   strike2 = col_double(),
##   ..   dip2 = col_double(),
##   ..   rake2 = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

Memahami isi dataset dan melihat jika ada missing value

summary(df)

##       tgl                  ot                lat               lon        
##  Min.   :2008-11-01   Length:89218      Min.   :-11.000   Min.   : 94.02  
##  1st Qu.:2015-09-20   Class1:hms        1st Qu.: -7.900   1st Qu.:113.65  
##  Median :2018-08-27   Class2:difftime   Median : -2.900   Median :121.19  
##  Mean   :2017-09-30   Mode  :numeric    Mean   : -3.396   Mean   :119.18  
##  3rd Qu.:2020-08-21                     3rd Qu.:  0.140   3rd Qu.:126.88  
##  Max.   :2022-10-31                     Max.   :  6.000   Max.   :142.00  
##                                                                           
##      depth             mag           remark             strike1     
##  Min.   :  2.00   Min.   :1.000   Length:89218       Min.   :  0.0  
##  1st Qu.: 10.00   1st Qu.:3.000   Class :character   1st Qu.:107.7  
##  Median : 16.00   Median :3.500   Mode  :character   Median :144.6  
##  Mean   : 49.32   Mean   :3.606                      Mean   :170.6  
##  3rd Qu.: 54.00   3rd Qu.:4.200                      3rd Qu.:219.6  
##  Max.   :750.00   Max.   :7.900                      Max.   :359.2  
##                                                      NA's   :86616  
##       dip1           rake1            strike2            dip2      
##  Min.   : 2.30   Min.   :-180.00   Min.   :  0.00   Min.   : 1.5   
##  1st Qu.:47.02   1st Qu.: -31.30   1st Qu.: 62.61   1st Qu.:39.3   
##  Median :62.30   Median :  57.60   Median :239.84   Median :58.5   
##  Mean   :60.22   Mean   :  30.14   Mean   :196.67   Mean   :56.6   
##  3rd Qu.:76.40   3rd Qu.: 100.60   3rd Qu.:296.92   3rd Qu.:74.7   
##  Max.   :90.00   Max.   : 180.00   Max.   :359.98   Max.   :90.0   
##  NA's   :86616   NA's   :86616     NA's   :86616    NA's   :86616  
##      rake2        
##  Min.   :-180.00  
##  1st Qu.: -20.58  
##  Median :  55.50  
##  Mean   :  34.44  
##  3rd Qu.: 112.00  
##  Max.   : 180.00  
##  NA's   :86616

lapply(df, class)

## $tgl
## [1] "Date"
## 
## $ot
## [1] "hms"      "difftime"
## 
## $lat
## [1] "numeric"
## 
## $lon
## [1] "numeric"
## 
## $depth
## [1] "numeric"
## 
## $mag
## [1] "numeric"
## 
## $remark
## [1] "character"
## 
## $strike1
## [1] "numeric"
## 
## $dip1
## [1] "numeric"
## 
## $rake1
## [1] "numeric"
## 
## $strike2
## [1] "numeric"
## 
## $dip2
## [1] "numeric"
## 
## $rake2
## [1] "numeric"

sum(complete.cases(df)) # checking on missing values

## [1] 2602

sum(complete.cases(df))/nrow(df)

## [1] 0.02916452

Membersihkan dan memanipulasi data

Merubah beberapa nama kolom

df2 <- df %>%
  rename(
    dates = tgl,
    times = ot,
    loc = remark
  )

Menambahkan kolom ‘months’ dan ‘years’ kedalam dataset

df3 <- df2 %>%
  mutate(
    years = year(dates),
    months= case_when(
                                              month(df2$dates)==01 ~ "January",
                                              month(df2$dates)==02 ~ "February",
                                              month(df2$dates)==03 ~ "March",
                                              month(df2$dates)==04 ~ "April",
                                              month(df2$dates)==05 ~ "May",
                                              month(df2$dates)==06 ~ "June",
                                              month(df2$dates)==07 ~ "July",
                                              month(df2$dates)==08 ~ "August",
                                              month(df2$dates)==09 ~ "September",
                                              month(df2$dates)==10 ~ "October",
                                              month(df2$dates)==11 ~ "November",
                                              month(df2$dates)==12 ~ "December"
                                            ))
colnames(df3)

##  [1] "dates"   "times"   "lat"     "lon"     "depth"   "mag"     "loc"    
##  [8] "strike1" "dip1"    "rake1"   "strike2" "dip2"    "rake2"   "years"  
## [15] "months"

Menghapus kolom yang tidak digunakan, serta melihat apakah masih ada data yang kurang lengkap

df4 <- select(df3, -c(strike1, dip1, rake1, strike2, dip2, rake2))
sum(is.na(df4)) # checking if still there are missing datas

## [1] 0

summary(is.na(df4))

##    dates           times            lat             lon         
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:89218     FALSE:89218     FALSE:89218     FALSE:89218    
##    depth            mag             loc            years        
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:89218     FALSE:89218     FALSE:89218     FALSE:89218    
##    months       
##  Mode :logical  
##  FALSE:89218

Melihat jika ada data duplikat

quakes_df <- df4 %>%
  distinct()

Analisa data

Plotting magnitudo dari sejumlah gempa bumi yang terjadi sepanjang tahun dari 2008 hingga 2022 (Oktober)

qplot(mag, data = quakes_df, bins = 100)

## Warning: `qplot()` was deprecated in ggplot2 3.4.0.

Mencoba melihat gempa bumi yang terjadi sejak 2008, dan terapkan dalam diagram batang berikut

year_df <- quakes_df %>%
  select(years) %>%
  group_by(years) %>%
  count() %>%
  arrange()
head(year_df)

## # A tibble: 6 × 2
## # Groups:   years [6]
##   years     n
##   <dbl> <int>
## 1  2008   316
## 2  2009  3534
## 3  2010  3437
## 4  2011  3940
## 5  2012  2422
## 6  2013  2719

# apply into bar chart
ggplot(data = quakes_df) +
  geom_bar(mapping = aes(x=years))+
  labs(title = "Earthquakes Incident Since 2008")

Apakah gempa terjadi berdasarkan bulan tertentu? skrip di bawah ini menunjukkan hasil selama bertahun-tahun sejak 2008

month_df <- quakes_df %>%
  select(months) %>%
  group_by(months) %>%
  count() %>%
  arrange()
head(month_df)

## # A tibble: 6 × 2
## # Groups:   months [6]
##   months       n
##   <chr>    <int>
## 1 April     7898
## 2 August    8451
## 3 December  7606
## 4 February  7304
## 5 January   7316
## 6 July      6895

# apply into bar chart
ggplot(data = quakes_df) +
  geom_bar(mapping = aes(x=months, fill=months)) +
  labs(title = "Earthquakes Incident by Month Since 2008")

Tentukan analisis per musim untuk melihat pola pada setiap musim pad masing-masing tahun

seasEq <- quakes_df %>%
  mutate(
    seasons = case_when(
      month(quakes_df$dates) == 10|
        month(quakes_df$dates) == 11|
        month(quakes_df$dates) == 12|
        month(quakes_df$dates) == 01|
        month(quakes_df$dates) == 02|
        month(quakes_df$dates) == 03 ~ "rainy season",
      month(quakes_df$dates) == 04|
        month(quakes_df$dates) == 05|
        month(quakes_df$dates) == 06|
        month(quakes_df$dates) == 07|
        month(quakes_df$dates) == 08|
        month(quakes_df$dates) == 09 ~ "dry season"
      
    )
  ) 

# show the scraped data
season_df <- seasEq %>%
  select(seasons) %>%
  group_by(seasons) %>%
  count() %>%
  arrange()
head(season_df)

## # A tibble: 2 × 2
## # Groups:   seasons [2]
##   seasons          n
##   <chr>        <int>
## 1 dry season   42526
## 2 rainy season 44090

# apply it into bar chart
ggplot(data = seasEq) +
  geom_bar(mapping = aes(x=seasons, fill=seasons)) +
  facet_wrap(~years) +
  labs(title = "Earthquakes Incident by Seasons")

Mencari pola jika kita mengkategorikan tingkat magnitudo ke dalam dataframe. Buat kolom baru bernama ‘mag_category’ yang membagi nilai di kolom ‘mag’ menjadi bagian dan menetapkan setiap nilai ke kategori yang sesuai (minor = 0-3,9 SR), ( sedang = 4-5 SR), (mayor = 5.1-7.9 SR), kemudian identifikasi pola yang ditunjukkan pada masing-masing tahun.

magEq <- quakes_df %>%
  mutate(mag_category = cut(
    mag,
    breaks = c(0, 3.9, 4.9, 5, 7.9),
    labels = c("minor", "moderate", "strong", "major")
  ))
head(magEq)

## # A tibble: 6 × 10
##   dates      times      lat   lon depth   mag loc           years months mag_c…¹
##   <date>     <time>   <dbl> <dbl> <dbl> <dbl> <chr>         <dbl> <chr>  <fct>  
## 1 2008-11-01 21:02:43 -9.18  119.    10   4.9 Sumba Region…  2008 Novem… modera…
## 2 2008-11-01 20:58:50 -6.55  130.    10   4.6 Banda Sea      2008 Novem… modera…
## 3 2008-11-01 17:43:12 -7.01  107.   121   3.7 Java - Indon…  2008 Novem… minor  
## 4 2008-11-01 16:24:14 -3.3   128.    10   3.2 Seram - Indo…  2008 Novem… minor  
## 5 2008-11-01 16:20:37 -6.41  130.    70   4.3 Banda Sea      2008 Novem… modera…
## 6 2008-11-01 14:47:00 -7.37  105.    18   3.3 Java - Indon…  2008 Novem… minor  
## # … with abbreviated variable name ¹mag_category

# plotting the finding (shown visual from data since 2008)
ggplot(magEq, aes(x = mag_category, fill = mag_category)) +
  geom_bar() +
  labs(title = "Number of earthquakes by magnitude category",
       x = "Magnitude Category",
       y = "Number of earthquakes") +
  facet_wrap(~years)

Mari kita coba filter data dengan magnitudo lebih dari 5 skala Richter yang menandakan bahwa gempa dengan kekuatan tersebut menimbulkan kerusakan dari ringan hingga parah, bahkan menimbulkan korban jiwa. Pertama, lihat apakah yang akan kita temukan jika menggunakan fungsi leaflet.

mag_maj <- magEq %>%
  filter(mag_category == "major")
leaflet() %>%
  addTiles() %>%
  addCircles(data = mag_maj,
             radius = mag_maj$mag,
             color = "red",
             fillOpacity = 0.1,
             weight = 1)

## Assuming "lon" and "lat" are longitude and latitude, respectively

Mengidentifikasi jenis gempa bumi berdasarkan kedalaman (ditampilkan data tahun 2008). Buat kolom baru bernama ‘eq_type’ yang membagi nilai di kolom ‘kedalaman’ ke dalam kotak dan menetapkan setiap nilai ke kategori yang sesuai (dangkal = 0-60 km), (menengah = 61-300 km), (dalam = > 300 km). Kemudian visualisasikan temuan tiap-tiap tahun.

depEq <- quakes_df %>%
  mutate(eq_type = cut(
    depth,
    breaks = c(0, 61, 300, 750),
    labels = c("shallow", "medium", "deep")
  ))
head(depEq)

## # A tibble: 6 × 10
##   dates      times      lat   lon depth   mag loc           years months eq_type
##   <date>     <time>   <dbl> <dbl> <dbl> <dbl> <chr>         <dbl> <chr>  <fct>  
## 1 2008-11-01 21:02:43 -9.18  119.    10   4.9 Sumba Region…  2008 Novem… shallow
## 2 2008-11-01 20:58:50 -6.55  130.    10   4.6 Banda Sea      2008 Novem… shallow
## 3 2008-11-01 17:43:12 -7.01  107.   121   3.7 Java - Indon…  2008 Novem… medium 
## 4 2008-11-01 16:24:14 -3.3   128.    10   3.2 Seram - Indo…  2008 Novem… shallow
## 5 2008-11-01 16:20:37 -6.41  130.    70   4.3 Banda Sea      2008 Novem… medium 
## 6 2008-11-01 14:47:00 -7.37  105.    18   3.3 Java - Indon…  2008 Novem… shallow

# visualize the finding
ggplot(depEq, aes(x = eq_type, fill = eq_type)) +
  geom_bar() +
  labs(title = "Number of earthquakes by Depth",
       x = "Earthquake Type",
       y = "Number of earthquakes") +
  facet_wrap(~years)

Mari kita lihat jika kita memetakan temuan menggunakan fungsi leaflet. Menentukan tipe kedalam ‘shallow’, karena berdasarkan riset gempa dengan kedalam dangkal paling merusak, menimbulkan bencana lain bahkan memakan korban jiwa.

depShl <- depEq %>%
  filter(eq_type == "shallow")
leaflet() %>%
  addTiles() %>%
  addCircles(data = depShl,
             radius = depShl$depth,
             color = "red",
             fillOpacity = 0.1,
             weight = 1)

## Assuming "lon" and "lat" are longitude and latitude, respectively

Meluangkan waktu sejenak untuk mengetahui korelasi apa antara kedalaman dan besarnya dan untuk memeriksa pengaruh antara kedua variabel ini, fungsi ‘stat_cor()’ digunakan untuk menghitung koefisien korelasi dan menambahkannya sebagai label teks ke plot.

ggplot(quakes_df, aes(x = depth, y = mag)) +
  geom_point() +
  stat_cor()

Mengelompokkan data berdasarkan ‘lat’ dan ‘long’ dan rangkum jumlah gempa bumi di setiap lokasi.

quake_locations <- quakes_df %>%
  group_by(lat, lon) %>%
  summarize(n = n())

## `summarise()` has grouped output by 'lat'. You can override using the `.groups`
## argument.

head(quake_locations)

## # A tibble: 6 × 3
## # Groups:   lat [1]
##     lat   lon     n
##   <dbl> <dbl> <int>
## 1   -11  114.     1
## 2   -11  115.     1
## 3   -11  116.     1
## 4   -11  117.     1
## 5   -11  118.     1
## 6   -11  119.     1

# Plot the data using a scatter plot
ggplot(quakes_df, aes(x = lat, y = lon)) +
  geom_point() +
  labs(title = "Number of earthquakes by location",
       x = "Latitude",
       y = "Longitude",
       size = "Number of earthquakes") +
  facet_wrap(~years)

# Create a map of the region
map_df <- filter(quakes_df, dates >= as.Date("2017-01-01") & dates <= as.Date("2022-10-31"))
leaflet() %>%
  addTiles() %>%
  addCircles(data = map_df,
             radius = map_df$mag,
             color = "blue",
             fillOpacity = 0.1,
             weight = 1)

## Assuming "lon" and "lat" are longitude and latitude, respectively

Temuan-temuan Penting

Setelah melakukan serangkaian proses analisa, ditemukan beberapa temuan-temuan yang dapat digunakan sebagai rencana mitigasi oleh pemangku kepentingan.

Kasus atau bencana gempa bumi, baik gempa kecil maupun besar terjadi pada tahun 2018, diikuti tahun 2019 dan tahun 2021.
Data juga menunjukkan sepanjang tahun 2008 hingga 2022, pada musim penghujan sering terjadi gempa bumi.
Sepanjang tahun 2008 hingga 2022, di Indonesia sering terjadi gempa minor atau gempa bermagnitudo anatara 0-3.9 SR, diikuti gempa sedang.
Banyak benacan gema bumi yang terjadi, baik gemap kecil maupun besar, dikategorikan sebagai gempa dangkal karena data mencatat kedalaman titik gempa tersebut berada dikedalaman antara 0-60 KM.
Data mencatat serta dikuatkan dengan hasil analisa menunjukkan bahwa lokasi yang sering terjadi gempa terjadi disepanjang pesisir barat Pulau Sumatera, sepanjang sisi selatan Pulau Jawa, Bali, Nusa Tenggara, wilayah Sulawesi dan kepulauan Maluku, serta sebagian Pulau Papua. Hal ini juga diimbangi dengan fakta bahwa daerah/kawasan-kawasan tersebut merupakan pertemuan 3 lempeng tektonik, yakni Lempeng Indo-Australia, Lempeng Eurasia dan Lempeng Pasifik.

Kesimpulan

Dengan melihat fenomena yang terjadi dapat disimpulkan bahwa sangat tidak mungkin bagi Indonesia untuk bebas dari bencana gempa bumi. Karena kondisi secara geografis, Indonesia berada pertmuan 3 lempeng tektonik yang mana setiap waktu terdapat aktivitas meskipun sekian milimeter. Hal yang bisa dilakukan adalah BMKG selaku pemangku kepentingan yang memantau adanya bencana, mensosialisasikan kepada masyarakat terutama yang tinggal dilokasi rawan gempa, berupa upaya penyelamatan jika terjadi bencana. Peringatan dini jika terjadi aktivitas dasar bumi yang bisa diakses masyarakat setiap waktu sebagai langkah preventif.
Selain itu, sebagai langkah preventif selanjutnya, pemerintah melalui sub-organisasi yang mengurusi tentang pembangunan bangunan, mensosialisasikan kepada masayarkat untuk mempertimbangkan menggunakan bangunan tahan gempa, terutama masyarakat yang tinggal dilokasi rawan gempa, seperti yang telah saya sebutkan diatas.

Analisa Gempa Indonesia 2008 - 2022

Nanda Akhirdianto

12/12/2022